Dependence of response functions and orbital functionals on occupation numbers 
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Explicitly orbital-dependent approximations to the exchange-correlation energy functional of den- 
sity functional theory typically not only depend on the single-particle Kohn-Sham orbitals but also 
on their occupation numbers in the ground state Slater determinant. The variational calculation 
of the corresponding exchange-correlation potentials with the Optimized Effective Potential (OEP) 
method therefore also requires a variation of the occupation numbers with respect to a variation 
in the effective single-particle potential, which is usually not taken into account. Here it is shown 
under which circumstances this procedure is justified. 



I. INTRODUCTION 

The central quantity of density functional theory p], Q , 
the exchange-correlation energy E xc , is a unique (though 
unknown) functional of the electron density. Popular 
approximations such as the local density approximation 
(LDA) and generalized gradient approximations (GGA's) 
express E xc as an explicit functional of the density. 

Recently, another class of approximations has at- 
tracted increasing interest: implicit density functionals, 
expressing E xc as explicit functionals of the Kohn-Sham 
single particle orbitals and energies and therefore only 
as implicit functionals of the density [1, 0] • Members of 
this class of functionals are the exact exchange functional 
(EXX), the popular hybrid functionals which mix GGA 
exchange with a fraction of exact exchange [EH Sill, the 
Perdew-Zunger self-interaction correction [9( and meta- 
GGA functionals [l(| HH, [HI which include the orbital 
kinetic energy density as a key ingredient. 

At zero temperature, the orbital functionals mentioned 
above depend on the occupied orbitals only. Other func- 
tionals, such as the second-order correlation energy of 
Gorling-Levy perturbation theory (l3j . in addition de- 
pend explicitly on the unoccupied orbitals and the orbital 
energies. Moreover, all these orbital functionals are not 
only explicit functionals of the orbitals but also explicit 
functionals of the occupation numbers which, in turn, 
depend on the single-particle orbital energies. This ad- 
ditional energy dependence is ignored in common imple- 
mentations of orbital- or energy-dependent functionals. 

In order to calculate the single-particle Kohn-Sham po- 
tential corresponding to a given orbital functional, the 
so-called Optimized Effective Potential (OEP) method is 
used [!, 0, [3] • The OEP method is a variational method 
which aims to find that local potential whose orbitals 
minimize the given total energy expression. In princi- 
ple, when performing the variation of the local potential 



one not only should vary the orbitals but also the orbital 
energies and occupation numbers. Typically, however, 
the variation with respect to the occupation numbers is 
not explicitly performed. In this work we will investigate 
when and why this is justified. 



II. DENSITY RESPONSE FUNCTION 

In this Section we analyze the problem of the eigen- 
value dependence of the occupation numbers in the den- 
sity and the non-interacting static linear density response 
function for various situations. We consider the case of 
zero temperature and distinguish between variations at 
fixed and variable particle number, i.e., for the canonical 
and grand-canonical ensemble. 



A. Fixed particle number 

The density of N non-interacting electrons (at zero 
temperature) moving in some electrostatic potential v s (r) 
is given by 



*(r)= j>i(r)| 2 , 



(1) 



where the single-particle orbitals are solutions of the 
Schrodinger equation 



(2) 



-— +v„(t) ) tpi(r) = enpi(r), 



and the sum in Eq. |T]) runs over the N occupied orbitals 
of the iV-electron Slater determinant. For the ground 
state density one can rewrite Eq. Jl} as 

n(r) = Y, 0(ep - £ 4 )k(r)| 2 = ]T fM(r)\ 2 , (3) 
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where the sum now runs over all orbitals. Sf is the 
Fermi energy, 9{x) is the Heavyside step function, and 
fi = 9(sF — Ei) is the occupation number of orbital fi(r). 
It is evident from Eq. ([3]) that the density not only de- 
pends on the (occupied) orbitals <ft(r) but also on the 
orbital energies £j, since the very specification of which 
orbitals are occupied and which unoccupied depends on 
their energies. 

Through Eq. ([2]), both of these quantities are func- 
tional of the potential v s (r), i.e., fi(r) — (pi[v s ](r), 
Ej = e,*[u s ]. The static density response function, which 
is the functional derivative of n with respect to v s , is 
therefore given as 



Sfi 



with 



The last step follows from first order perturbation theory, 
which can be used to obtain 



<V»(r) 
Sv s (r>) 



J2 l Pk( F ) l Pt( r ') ( Pi( r ') 



(6) 



For simplicity, we have assumed a non-degenerate single- 
particle spectrum. 

Usually, x( r i r ') of Eq. ((5]) is taken as the static density 
response function instead of \. Both expressions differ by 
the first term on the right hand side of Eq. ((4]) , becoming 
identical only if this term vanishes. In order to see when 
and how this happens we consider two cases. 

Case 1 comprises systems for which the single-particle 
spectrum has a finite gap between the highest occupied 
orbital (eigenvalue ejv) and the lowest unoccupied orbital 
(eigenvalue ejv+i)- Then the Fermi energy sf lies strictly 
between these two orbital energies, Sn < £f < £n+i- 
Within the single-particle gap, the position of sf is arbi- 
trary (at zero temperature) . The important point now is 
that upon (infinitesimal) variation of the potential v s , ep 
remains fixed and does not need to be varied. The rea- 
son is that the variation 5sn of £jv due to the variation 
of v s is infinitesimal as well and ef can be chosen such 
that £f > £n + $£n, thus leaving the particle number 
unchanged. Then the functional derivative of the occu- 
pation number with respect to v s becomes 



Sfi d9(e F - Ssi 



Sv s (r) 



de l Sv s (v) 



= -5{e F ~e l )\^{v)\\ (7) 



where 8(x) is the Dirac delta function and we used the 
relation 



Sv s (r) 



(8) 



which can be obtained from first-order perturbation the- 
ory. In the present case, the Fermi energy (which is in 
the single-particle gap) is not equal to any of the single- 
particle energies, the delta function in Eq. ([7]) vanishes 
and x(r, r ') °f Eq. ((4]) coincides with the usual form of 
the static density response function of Eq. ([5|) . 

Case 2 is the case of a vanishing single-particle gap, i.e., 
the case of an open-shell or metallic system. For nota- 
tional simplicity, in the following discussion we still work 
with the assumption of a non-degenerate single-particle 
spectrum. Of course, particularly for open-shell systems, 
this assumption is inappropriate. The more general case 
including degenerate single-particle orbitals is discussed 
in Appendix [Al 

The crucial difference to case 1 is that an infinitesi- 
mal variation of the potential v s now not only leads to 
a variation Sei of the single-particle energies but also to 
a variation Ssf of the Fermi energy. This latter varia- 
tion has to be taken into account in order for the particle 
number to be conserved (i.e., the infinitesimal variation 
SN of the particle number upon variation of the poten- 
tial strictly has to vanish, SN = 0). Then the functional 
derivative of the occupation number with respect to the 
potential consists of two terms and reads 



Sfi _ d9(e F - £i) Se F | d9(e F - £%) Sej 
Sv s (r) de F Sv s (r) 

= <5(£ F -£,)(l^(r)| 2 



I^WI 2 ) 



Sv s (r) 



(9) 



where <pp is the highest occupied orbital with orbital en- 
ergy equal to the Fermi energy. Due to the delta func- 
tion, the r.h.s. of Eq. |9]) vanishes and again x(r,r') of 
Eq. (|U) coincides with the static density response func- 
tion %(r,r') of the form given in Eq. (J5]). 

From Eq. (f4]) the linear change in the density due to the 
perturbation Sv s (r) is Sn(r) = fd 3 r'x(r,r')5v s (r'). One 
can then check explicitly that the result x(r, r') = x(r, r') 
obtained here is fully consistent with a fixed number of 
particles: 



SN = Jd 3 r5n{r)= jd 3 r'Sv s (r') y d 3 r x(r, r') 

= jd 3 r'Sv s (r') Jd 3 r x (r,r') = 0, (10) 

where the last equality follows from the orthonormality 
of the single-particle orbitals. 

B. Grand canonical ensemble 

The analysis is slightly altered if the system of non- 
interacting electrons is connected to a particle bath, 
i.e., for the grand canonical ensemble characterized by 
a chemical potential \x. The density (at zero tempera- 
ture) is then given by 



l (r) = ^^M- £ i)ki(r)| 2 = ^/ i |^(r)| I 



(11) 



3 



where the occupation number now is given by /, = 
— £j) and the sum again runs over all single-particle 
orbitals. When varying the occupation numbers with re- 
spect to variations of the potential, the chemical potential 
remains constant, independent of the single-particle spec- 
trum having a finite or vanishing gap at /x. The variation 
of fi then is obtained similarly to case 1 of the previous 
subsection as 



Sfi 



d6(fx - Ssi 



Sv s (r) 



ds;, 



Sv s (r) 



-£( M - £l )l^(r)| 2 . (12) 



This term does not vanish if the chemical potential is 
aligned with one of the single-particle energies and the 
static density response function for the grand-canonical 
ensemble reads 



X(r,r') = x(r,r') e l )\(p l {r){' 



(13) 



It is worth noting that now, due to the second term 
on the r.h.s. of Eq. (13]), SN (Eq. JTUJ) ) is different from 
zero which is of course consistent with the fact that here 
we are dealing with an open system. 



III. IMPLICATIONS FOR THE OPTIMIZED 
EFFECTIVE POTENTIAL 



The central idea of density functional theory is to write 
the ground state energy E to t of N interacting electrons 
moving in an external electrostatic potential Vq(t) as a 
functional of the ground-state density. This energy func- 
tional may then be split into various pieces as 

Etot = T s [n] + Jd 3 r v {r)n(r) + U[n] + E xc [n], (14) 

where T s [n] is the kinetic energy functional of non- 
interacting electrons, 



U[n] = \ Jd\ JS 



,ra(r)n(r') 



(15) 



is the classical electrostatic energy and E xc is the 
exchange-correlation energy functional which incorpo- 
rates all complicated many-body effects and in practice 
has to be approximated. Minimization of Eq. (fl"4"|) with 
respect to the density leads to an effective single-particle 
equation of the form of Eq. where the effective po- 
tential is 

v s (r) = «o(r) + jd\'^- + v xc {v), (16) 
J |r-r'| 

with the exchange-correlation potential 

SE XC 



v xc (r) 



Sn(r) 



(17) 



While the most popular approximations to the 
exchange-correlation energy E xc are explicit functionals 



of the density, there has been increasing interest in an- 
other class of approximations which are are only implicit 
functionals of the density. These functionals instead de- 
pend explicitly on the Kohn-Sham single-particle orbitals 
as well as on the Kohn-Sham orbital energies. One ex- 
ample for such a functional is the exact exchange energy 
given as 



E 



EXX 



where 



7(r,rO = y>¥>i(rK(r / ) 



dV 



,l7(r,r')| 5 



(18) 



(19) 



is the single-particle density matrix. As one can see, 
E x xx depends on the single-particle energies through 
the occupation numbers /j. Other functionals such as, 
e.g., the correlation energy functional of second-order 
Gorling-Levy perturbation theory [l3[ , depend on the or- 
bital energies also in other ways (see below). 

In order to distinguish a genuine dependence on orbital 
energies from a dependence on occupation numbers we 
write for a general exchange-correlation energy functional 
E xc = E xc [{(pi},{ei},{fi}]. The exchange-correlation 
potential of such a functional can be computed by us- 
ing the chain rule of functional differentiation as 



v xc (r) 



SE XC 
Sn(r) 



, , SE Tr SvJr') 
dV *,.£) 4) ' (20) 



Acting with the density response operator on both 
sides of this equation one arrives at 



dV^ c (r')x(r',r) 



d 3 r ' 5E *c 



6v s (r>) 




5E XC 



Of, 



5<pi(r>) 
Sfi 



<W) 

Sv s (r) 



+ c.c. 



(21) 



In the last step we have used the chain rule once again 
and we also emphasize in the notation that when varying 
with respect to one set of variables (orbitals, orbital en- 
ergies or occupation numbers) the other variables remain 
fixed. 

Eq. (|2Tj) is the OEP integral equation in its general 
form. For a given approximate E xc , this equation defines 
the corresponding v xc (r) and has to be solved in a self- 
consistent way together with the Kohn-Sham equations 
(Eq. (HI)). It differs in three ways from the form most 
commonly found in the literature (see, e.g., Refs.H,0and 
references therein). One, the explicit energy dependence, 
is handled in a similar way as is the orbital dependence, 
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via the chain rule. The other two arise from the implicit 
energy dependence of the occupation numbers, and are 
our main concern here. Similar to the discussion in the 
previous section we will again distinguish between the 
two cases of fixed particle number and systems in contact 
with a particle bath and discuss the role of these extra 
terms in both cases. 



A. Fixed particle number 

As we have seen in Section [Til f° r the case of fixed par- 
ticle number at zero temperature the functional deriva- 
tive 6fi/8v s (r) vanishes both for systems with a finite 
and vanishing HOMO-LUMO gap. This has two conse- 
quences for Eq. (|2ip : first, we can replace the response 
function x by the function x of Eq. ([5]) and second, the 
last term on the r.h.s. of Eq. ([2"Tj) drops out. Therefore, 
the OEP equation reads 



dV^ c (r')x(r',r) 




{e k },{M Svs ^ 



(22) 



This equation shows that despite the dependence of E xc 
on the occupation numbers (which, in turn, depend on 
the orbital energies), the variation with respect to these 
occupation numbers may be omitted for the calculation 
of the OEP integral equation for the exchange-correlation 
potential. This is, of course, what has been done in the 
vast majority of cases discussed in the literature. 

We note in passing that integrating Eq. (l2"2")) over all 
space and using the orthornormality of the Kohn-Sham 
orbitals one can deduce the sum rule 11511 



E 



9E XC 



de. 



= 



Wk},{fk} 



(23) 



On quite general grounds, one expects that for an iso- 
lated system with a fixed number of particles, v xc {r) is 
only defined up to a constant. To check if Eq. f2"2"j) meets 
this condition we need an explicit expression for E xc . As 
a non-trivial example, we use 

E xc n E* xx + Ej 2 \ (24) 

where E x xx is the exact exchange energy of Eq. (fT5)) and 

(2) 

Ec is the second-order correlation energy of Gorling- 
Levy perturbation theory [l3|, [l^, [TtJ defined by 

E c {2) = E Ctl + E ca , (25) 



where 



E r 



E 



- fj) 

(Si - Sj) 



l(*kb'} + E/^*')l 2 ' ( 26 ) 



and 

E = ly M(l-A)(l-/z) 
c ' 2 2 Z-f (e 4 + Sj -e k - ei) 

(ij\\kl)[(kl\\ij)-(kl\\ji)}. (27) 
In the equations above we have used the notations 



mkl)= /d 3 r /dV wl'^Ww-cowM 



r — r 



and 



(i\v x \j) = ld 3 r ipi*(r)v x (r)(pj(r) 



(28) 



(29) 



Suppose now that we introduce a rigid shift t! s (r) — ► 
v s (r) + C in the effective single particle potential of 
Eq. As a result, if {<fi}, {£%}, {fi} are a set of solu- 
tions for V s (r), the solutions for v s (r) + C are {(fii} 7 {ei + 
C},{fi}. This holds provided that Eq. (f2"2"]l determines 
^xc(r) only up to a constant. Inspection of Eq. (j22j) con- 
firms that this is the case: the l.h.s. is invariant under a 
rigid shift of v xc (r), and Eqs. (|18[) and ([23]) are invariant 
under the change {ei} — * {ei + C}. 



B. Grand canonical ensemble 

The situation is different if the system is in con- 
tact with a particle bath. Since in this case 6fi/6v s (r) 
does not vanish one has to use the full OEP equation 
(f2"Tj) . Here the dependence of both the density and the 
exchange-correlation energy on the occupation numbers 
has been taken into account explicitly when performing 
the variations and the two extra terms resulting from this 
variation cannot be neglected. Applications of this OEP 
formalism for open systems have been reported for quasi 
two-dimensional electron gases (2DEG) in n-doped semi- 
conductor quantum wells where the n-doped regions act 

as particle reservoirs mis Eg. 

As another consequence of the extra terms, integration 
of Eq. I|2ip over all space leads to the modified sum rule 



J^<5(^ - Ei)v XC:l = ^ 

i i 

dE x 



dE x 



de. 



5{pL - Ei) 



where 



{Vfc},{£fc} 



d 3 rv xc (r) \(pi(r)\ 2 



(30) 



(31) 



We take the exact exchange functional (fT5)) as an example 
for a functional which does not explicitly depend on the 
single-particle energies. In this case, the first term on 
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the r.h.s. of Eq. (|30|) vanishes. If there exists a single- 
particle state whose energy equals the chemical potential, 
En = Mi we then obtain 



-EXX 

J x,N 



dE. 



EXX 



df. 



N 



(32) 



This relation is the complete analogue for the grand 
canonical ensemble of a well-known relation for fixed par- 
ticle number which reads [2l], [22|, HH 



-EXX 
'x.N 



-EXX 
l x.N 



where 



-EXX 



= ~T fd 3 r^ N (r) 

JN J 



6E. 



EXX 



5(p N (r) 



(33) 



(34) 



For open 2DEG's, this relation has been obtained previ- 
ously by studying the asymptotic behavior of the exact- 
exchange potential [l9(. 

For the grand-canonical ensemble, v xc (r) is fully de- 
termined by Eq. (|2ip since this equation is not invariant 
under a rigid shift of the potential: the l.h.s. is not in- 
variant due to the extra term in x(r, r') in Eq. (fT3|) . The 
r.h.s. is not invariant because E xc changes under the 
transformation {e,;} — * {ei + C}. This is due to the fact 
that the chemical potential \x (which is determined by the 
particle reservoirs) remains fixed in the grand canonical 
ensemble and the above transformation leads to a change 
in the set of occupation numbers and self-consistent KS 
orbitals, {f{\ and {</?i}, respectively. 



IV. CONCLUSIONS 

In this work we have addressed the question why and 
when one can ignore the explicit dependence on the 
orbital occupation numbers (which in turn depend ex- 
plicitly on the orbital energies) when calculating both 
the static linear density response function and the effec- 
tive single-particle potential corresponding to an orbital- 
dependent exchange-correlation energy functional. Wc 
have shown that the variation of the occupation numbers 
may safely be neglected for systems with fixed particle 
number. For systems connected to a particle bath, how- 
ever, this variation leads to non-vanishing contributions 
and needs to be taken into account. 



APPENDIX A: DEGENERATE 
SINGLE-PARTICLE SPECTRUM 

In general, the single-particle spectrum will have eigen- 
values which may be degenerate. In particular, in the 
case of open-shell systems, the energy of the highest oc- 
cupied orbital is degenerate and the arguments of the 
Case 2 discussed in Section III Al need to be modified. 



As degeneracy is in almost all cases related to symme- 
try we will use the language of group theory. In partic- 
ular, the single-particle orbitals will be labelled by the 
a complete set of quantum numbers {n, I, m} where n is 
the principal quantum number (which is not related to 
symmetry), I is a label denoting the irreducible repre- 
sentation of the symmetry group Q of the single-particle 
potential V s (r), and m labels a partner within that rep- 
resentation. The single-particle equation now reads 



2 



+ v s (r) (p n im(r) = e n np n i m (v) 



(Al) 



and it should be noted that the eigenvalue e„i is inde- 
pendent of the partner label m. Furthermore, writing 
the energy eigenvalue as a functional of the potential, 
£nz[v s ], one has to keep in mind that this functional is 
only well defined for potentials which are invariant under 
the transformations of the symmetry group Q because 
I refers to an irreducible representation of that group. 
Therefore, we calculate the variation of the orbital en- 
ergies, Se n i = e n i[v s + Sv s ] - e n i[v s ] resulting from a 
variation Sv s (r) which preserves the symmetry of v s (r). 
Replacing v s -> v s + Sv s , ip n i m -> tp n i m + Sip n i m , and 
£ni — > £ni + 5e n i in Eq. (|A1|) . one finds that the first- 
order change in the energy eigenvalue is given by 



8e r 



& z r\ip nlm {v)\ 2 5v s {v) 



(A2) 



Summing this equation over the partner label m one ob- 
tains 

d„i5e nl = /d 3 r^ \tp rdm (r)\ 2 Sv s {r) , (A3) 



where d n i is the degeneracy of e n i ■ Now we note that the 
single-particle orbitals </? ra im(r) may be written as 



Vnim(r) = R n i(r)Xi m (r) , 



(A4) 



where R n i(r) is a totally symmetric function which is 
invariant under all symmetry transformations T of the 
group Q and A; m (r) is a function which transforms ac- 
cording to the irreducible representation / of Q, i.e., 



X lm (BT\T)T) =^r«(T) mIm X im ,(r) 



(A5) 



Here, R{T) is a 3 x 3 matrix describing the symme- 
try operation T S Q in three dimensional space and 
T"'(T) is the representation matrix of group element T 
in the irreducible representation I of Q. Noting now that 
J2 m l^m( r )| 2 1S a totally symmetric function, we find for 
the functional derivative 



Sv s (r) 
where we have defined 

R nl (r) = -4= 



\Rm(r)\ 2 , 



Rni(r) t 



\X lm (r) 



(A6) 



(A7) 
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which is again invariant under all symmetry transforma- 
tions of the group Q. 

Eq. (|A6|) will shortly be used to repeat the arguments 
of Section III Al for the degenerate, open-shell case. Before 
we do so, we point out that the definition of the density of 
Eq. ([3]) needs to be modified because not all orbitals with 
energy ep are (fully) occupied. This can be achieved, e.g., 
by writing the density as 



and 



n(r) 



fnlm\<finlm(r)\' 



(A8) 



and occupying all degenerate orbitals of the partially 
filled subshell with the same fractional number of elec- 
trons, i.e., by defining the occupation number of the par- 
tially filled subshell by f n i m = / ni = (n n i/d„i)9(£p - s„i) 
where n n i is the number of electrons in the open subshell. 
With this definition, the static density response function 
reads 



X(r,r') = x(r,r') + ^ 



Sfnl 

Sv s (r) 



Wnlr, 



(A9) 



with 



X(r,r') 



E /» 



n,l,rr, 



' fn'i'm' ( r )<Pn>i'm' (r')yw (r)<; m (r') +cc ^ ^ AlQ j 



Sfnl _ n n i 
Sv s (r) d n i 

(\R F (r)\ 2 



8(SF - Enl) 
-\Rnl{v)\ 2 



(All) 



where the last equality follows because the total symmet- 
ric part of degenerate orbitals is identical. Therefore, just 
as in the non-degenerate case at fixed particle number, 
the functional derivative w.r.t. the occupation numbers 
may be neglected both in the calculation of the density 
response function as well as in the derivation of the OEP 
equation. 
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